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Motivated by the study of rare events for a typical genetic switching model in systems biology, in this paper 
we aim to establish the general two-scale large deviations for chemical reaction systems. We build a formal 
approach to explicitly obtain the large deviation rate functionals for the considered two-scale processes based 
upon the second-quantization path integral technique. We get three important types of large deviation results 
when the underlying two times scales are in three different regimes. This is realized by singular perturbation 
analysis to the rate functionals obtained by path integral. We find that the three regimes possess the same 
deterministic mean-field limit but completely different chemical Langevin approximations. The obtained 
results are natural extensions of the classical large volume limit for chemical reactions. We also discuss its 
implication on the single-molecule Michaelis-Menten kinetics. Our framework and results can be applied to 
understand general multi-scale systems including diffusion processes. 

Keywords: two-scale large deviations, second quantization path integral, mean field limit, chemical Langevin 
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I. INTRODUCTION 

In recent years there has been a growing interest in studying the rare transitions for fast-slow stochastic dynamics in 
biology n n uni HU nil na [n 1311 and 133 ] • In computational neuroscience, the stochastic hybrid system is utilized to 
model the fast switching of ion channels and the membrane voltage evolves according to a dynamics which depends on 
the ion channel states. In systems biology, people are interested in the phenotypic switching of the cells modeled by the 
central dogma, which involves fast switching of DNA states between active and inactive states and the transcriptional 
and translational processes with different rates depending on the DNA states. In both cases, the transition rates and 
the most probable transition paths between different stable fixed points are issues being investigated in the literature. 
The main approaches include the WKB asymptot ics and the path integral formulations. However, mathematically 
it falls in the field of large deviation theory (LDTjPEsElHlSl and the rigorous results for these types of problems are 
very limitecP^. It is also meaningful t o remar k that there is a close connection between the LDT and the popular 
landscape theory for biological systemJi^IIlIlS] 

In this paper, we will continue our program to study the two-scale large deviations for c hemic al kinetic systems. To 
illustrate our points more concretely, let us consider a canonical genetic switching modell^*^ in systems biology as 
shown in Fig. Dynamics of this self-regulating genetic system can be described by the following chemical master 
equation 

dtP{n,t) = ^ (E-i - l)P(n,t) +7(E1 - l)[nP(n,t)] -f ^ P{n,t), ( 1 ) 

where the two-component vector P{n,t) = {Po{n,t), Pi{n,t))'^, and Pj{n,t) is the probability distribution function 
that the system has n protein copy numbers in the DNA active (j = 1 ) or inactive state (_) = 0 ). The raising or 
lowering operator Ej) is defined through E^/i(n) = h{n + k) for any function h depending on n. kg and ki are protein 
synthesis rates, 7 is the degradation rate constant, and /(«•), g{n) are switching rates between two DNA states. 

The biologically relevant parameter setup is 7 ~ 0 ( 1 ) and fco/7, ^1/7 both large. We will not consider more 
detailed regimes concerning the magnitudes of ki and ko although one usually has fci 3 > k^ in realistic situations. 
This does not affect the main point in this paper. In this case, the average number of proteins at steady state is of 
order fci/7. Now let us define the small parameter e « 7/^1 or 7/fco, thus the characteristic number of proteins is 
n ^ 0 (e“^). In our fast-slow genetic switching model, we define the switching rates f{n),g{n) ^ 0 (e““), and the 
realistic situations can be classified into the following three typical regimes: 
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FIG. 1: A typical fast-slow genetic switching model considered in systems biology. Left panel: Schematics of the 
chemical reaction schemes, where the switching rates f{n), g{n) are usually large. Right panel: Direct Monte Carlo 
simulations of the genetic switching model. The bi-stability is clearly observed from the time series of protein copy 

numbers. 


• Case 1: a > 1, i.e. the genetic switching process is much faster than the translation process; 

• Case 2: a = 1, i.e. the switching rates are comparable to the translation rates; 

• Case 3: 0 < a < 1, i.e. the translation process is much faster than the genetic switching process. 

In Case 2, the WKB asymptotics and the rigorous LDT results have been established for a similar model which 
takes into account the mRNA fluctuatiorfiS. The obtained LDT rate functional is utilized to find the most probable 
transition path and characterize the rate of transitions between the high and low expression states. Furthermore, the 
authors have shown that the Hamiltonian obtained from LDT is convex with respect to the momentum variable, which 
is one key point in designing robust numerical algorithms. In Case 3, the researchers typically take the continuum 
limit to the translation process at first since it is even faster than the switching procesa^^. With this approach, one 
obtains a stochastic hybrid system which resembles similar form as those for ion channels considered in computational 
neuroscience. So far, the WKB asymptotics and path integral formulations are both proposed for stochastic hybrid 
systems. The Case 1 is also studied with WKB asymptotics applied to the averaged system with respect to the fast 
switching process. 

From the authors’ point of view, the approaches employed in m are like taking a repeated limit to the switching 
and translation processes according to their relative magnitudes. More concretely, when DNA switching is much 
faster than the protein synthesis, the equilibrium pre-averaging of the switching process is taken in [10] at first and 
one gets a pure translation process with effective translation rates; however when the protein synthesis is much faster 
than DNA switching, the large volume limit is taken to the translation process at first and one gets a stochastic hybrid 
systenJi^. Similar ideas and techniques are adopted in [20] and [20] as well, which discussed different timescale issues 
for the gene expression model. With this understanding, it will be interesting to investigate the double limit of the 
original process instead of taking average with respect to one faster process at first. Mathematically it is also desirable 
to establish the large deviations for the original system with two time scales but different magnitudes. In fac t, it is 
the main motivation of this paper. We will utilize the Doi-Peliti second quantization path integral formalisrrPI^lISIl 
to study the general two-scale large deviations for the genetic switching models. As we will see, although the second 
quantization path integral for the spin-boson type model is formal, it is an effective approach to derive the large 
deviation results for chemical jump processes. Compared with the classical path integral formalism for diffusion 
processes, the second quantization path integral for chemical jump processes formulates the weight of each path in 
an extended space which involves both coordinate and momentum variables. This makes that the large deviation 
result can be given through a Hamiltonian with explicit formula, which resolves the dilemma that the Lagrangian 
in the rate functional does not have a closed form. This is important for further theoretical and numerical studies. 
Mathematically, rigorously establishing the LDT obtained from the formal approach in this paper is in progress based 
on our previous analysis m- 

Let us briefly illustrate our general two-scale LDT results. We will show that the Lagrangian obtained from the 
second quantization path integral comprises of two parts, which correspond to the switching and translation processes, 
respectively. However, what we are interested in is the LDT only for the concentration of proteins. The different 
magnitudes of the switching and translation rates essentially lead to a singularly perturbed variational problem, 
which has different dominant terms and different scaling limits in the cases of 0 < a < 1 and a > 1. When a = 1, 
the Lag rangians from both parts contribute equally, and we get a result which combines the Donsker-Varadhan type 
LDli^for the occupation measure of DNA states and the large volume type LDT^^for the small noise perturbation 
altogether m- As the LDT gives the sharpest characterization of the considered two-scale chemical kinetic system, 
we can obtain the deterministic mean field ODEs and the chemical Langevin approximation for the system based on 
the local analysis of the large deviation result^. This corresponds to the law of large numbers (LLN) and the central 
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limit theorem (CLT) for the process. We found that the three cases possess the same mean field ODEs. However, the 
chemical Langevin approximations for them are quite different. If a > 1, only the fluctuation from protein translation 
process survives. If 0 < a < 1, only the fluctuation from genetic switching process survives. And if a = 1, both 
fluctuations from protein translation and renetic switching processes contribute. Similar results are also valid for the 
single-molecule Michaelis-Menten kineticJ^l -^^ith slight modifications (c.f. Section]^. Our study extends the insights 
about the chemical kinetic systems in the classical large volume limit, and the methodology we introduced here can 
be applied to other multiscale problems in many fields. 

The rest of the paper is organized as follows. In Section [TTj we introduce the chemical master equation and apply 
the Doi-Peliti path integral formalism to the considered model. We then rescale the system with system size e~^ 
and get the abstract LDT result based on singular perturbation analysis in Section m In Section [TV] we apply our 
abstract result to the two-state genetic switching model and present the mean field limits and chemical Langevin 
approximations. In Section|Vj we apply our result to the well-known single-molecule Michaelis-Menten kinetics and 


mention the implications. Finally we make the conclusion and related discussions in Section VI 


II. TRANSITION PROBABILITY IN A PATH INTEGRAL FORM 


We start from a more general model rather than Eq. Q. Assume that the DNA switching could occur among N 
possible states (N = 2 for the model shown in Fig. and the chemical master equation (CME) for the biological 
reaction network reads 




P(n,t) 


P(n — l,t) — P(n, t) 

(n + l)P{n -I- 1, t) - nP{n, t) 
-|-Q'^P(n, t). 


( 2 ) 


Here P{n,t) = (Pi(n,t ),..., PNin,t))'^, Pj{n,t) is the probability distribution function that the system has n protein 
copy numbers and the switch is in state j at time t. A is a diagonal matrix with diagonal entry kj as the protein 
synthesis rate in state j. 7 is the protein degradation rate and Q = is the transition rate matrix among 

different DNA states. Thus qjk{n) > 0 for any j k and Qjk{n) = 0. We assume that the switching process is 

ergodic. _ 

Now we follow the Doi-Peliti’s approach to establish the path integral formalism of the CME (j^ 7 | 24 | 3i |32] 
the creation, annihilation operators a\a and the state function \^p) as 


J. % ^ 00 

a' \n) = |n -f 1), a |n) = n |n — I) and \^f)) = } P{n, t) \n). 

Then the CME ([^ can be written in a second-quantized form 

dt lip) = fl\^p), (3) 


where the operator 


r2 = A{a' — 1) -I- 7(0 — a^a) -f (4) 

and Q is obtained from Q by replacing the transition rates qij{n) with operators qijia^a). 

From Eq. ([^, the transition probability P(n/,r|ni,0) of finding product copy number Uf at time t = t starting 
from rii at t = 0 has the form 


P{nf,T\ni,0) = {nf\exp{flT) \n,) 

= {iT'fl exp(r2At)’’/^‘ jn^), (5) 

where \ni) = {\ni), \ni) , ■ ■ ■ ,\ni))’^ is an TV-dimensional column vector and (n/| = ((n/j,-- - ,{nf\) is an N- 
dimens ional row vector. Following Zhang et al.123, we utilize the coherent state representation and a resolution of 
identitji^i^ as 
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where 


with 


\^P^) = 


^ros2^p*‘>i/2 

N^ob 


2- r*n«;^ ^ p'^4>N /2 




|z) = e“ " |0) , (z| = (0| e“^ z = ne"*^, z = e*^, 


and i is the imaginary unit. The variable n has the interpretation that it characterizes the mean protein number 
in the coherent states. Define Cj = {2/N)cos^ {9j/N). The Cj gives the occupation probability of DNA at state j 
from the probabilistic interpretation of quantum mechanics. They satisfy the normalization condition ~ 

Correspondingly the phase variable (p can be chosen to satisfy 4'j = 0 for convenience. These choices are 

consistent with the resolution of identity. As a consequence, we have only iV— 1 independent unknowns (ci,..., cn-i) 
in c = (ci,..., Cat), which is equivalent to use 9 = (0i,..., 07 v), and A^ — 1 unknowns [pi,... ,4 >n-i) in (/> = 

(^ 1 , . . . ,(j)N)- 

Inserting § into (§, the transition probability density can be represented as a path integral form 

DnDjjDcDcpexp (^~ J Ldt^ , ( 7 ) 


P(n/, T|ni, 0 ) = Const. X J 


where the Lagrangian L is defined as 


. N-l 
. CLTl . 

L = il3— + ^ 


dt 


i=i 


djcN - Cj)/2 

' dt 


- H{n,c,i/3,icp), 


( 8 ) 


and the Hamiltonian 


H{n, c, if3, i(f)) = Hi{n, c, i/3) + H 2 {n, c, icp) 


(9) 


with Hi for the translation process 


N 

Hi{n,c,i/3) = kjCj [exp(f/3) — 1] + 'yn[exp(—i/3) — 1] 
j=i 


( 10 ) 


and H 2 for the switching process 


N 


H2{n,c,i(f)) = - 1 ). 


( 11 ) 


In Eq. 0. the outer 4-fold integral is taken in the path space with respect to n{t), I3(t), c{t) and (f>{t), which are full 
trajectories in [0,r]. The terms involving \ni) and (n/| have been absorbed to the constant before the integral. The 
path integral formulation Q makes the weight of each trajectory explicit. The form of Lagrangian (|^ suggests the 
interpretation that the pairs i/3 and n, ip and c are conjugate variables. 

To study the associated LDT, we must have a small parameter e and a deterministic limit as e —^ 0. This could be 
chosen as the inverse of typical system size e = y/Zci = V~^. As stated in the introduction section, we assume 


kj 


Qij 


a > 0 


and define 

X = ne, kj = kje, qjk{x) = qjk{n)t 
With these definitions, Eq. 0 can be rewritten as 


P{nf, rim, 0) = Const, x 


DxD/3DcDc/)exp ( — 


Lidt - 


Lodi 


( 12 ) 


(13) 


(14) 
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where the rescaled Lagrangian 


and rescaled Hamiltonian 


~ dx ~ 

Li = 7 / 3 — - Hi{x,c,il3), 


N-l 

L2 = e°‘Yl 

i=i 


d{cN -Cj)/2 

dt 




N 

Hi{x,c,ij3) = fcjCj[exp(z/3) — 1] + 7 a:[exp(—z/3) — 1], 
i=i 

N 

H2{x,c,i4>) = - 1 ). 


(15) 

(16) 


(17) 

(18) 


Using the method of steepest descent asymptotics, the integration over /3 and 4> can be approximated by simply using 
the value of the integrand at the saddle pointP. Thus, we get 


P{nf,T\ni,0) 


DxDc exp 



Li{x, X, c)dt 


^ L2{x,c)dt^ 


(19) 


where 


Li(x, X, c) = supjpi; — Hi{x, c,p)}, 

p 

L2{x, c) = sup{-H2{x, c, v?)}. 


( 20 ) 

( 21 ) 


Note the term tp ■ c does not appear in Eg . (|21[ ) because of the factor e“ in the first term of Eg. (16). Formally, the 
functional appearing in the exponential in ( |19[ ) is a competition between the rate functional f Lidt which corresponding 
to the translation process and the rate functional f L 2 dt which corresponding to the switching process. It is interesting 
to observe that the Lagrangian Li corresponds to the large volume type LDT rate function for the small noise 
perturbatio iP^ a nd L 2 corresponds to the Donsker-Varadhan type LDT rate function for the occupation measure of 
DNA state^^mH, xhe second guantization path integral perfectly reveals the intrinsic structure of the considered 
two-scale chemical kinetic process. 


III. FORMULATION OF THE LDT IN A GENERAL SETTING 


The transition probability (191 contains the LDT information about the variables x and c. However in most cases, 
one is only interested in slow variables, i.e. the concentration of protein in our case, which is also the observable in 
experiments. In this sense, we must integrate over c-space. It turns out the final result depends on the value of a and 
we will have three typical regimes. In what follows, we will discuss different outcomes in different regimes separately. 

(i). Case 1: a > 1. The switching process is much faster than the translation process. 


In this case, we can rewrite Eg. (191 as 


P(n/, r|ni, 0) oc J DxDcexp — - J dt ^Li{x,x,c) + ^_^ L 2 {x,c)^'^ . 


( 22 ) 


To integrate over c-space, we take the Laplace asymptotics for each t. The Lagrangian for x has the form 

Lx{x,x) = inf ■|Li(a:,i:,c) -|- ^_^ L 2 {x,c)'^ as e —)■ 0-I-. (23) 

Since L 2 > 0 and » 1, the term e^~°'L 2 {x,c) dominates. From the assumption that the switching process is 
ergodic, for a given x, L 2 {x,c) achie ves its minimum L 2 {x,c) = 0 at a single point c = cq(x), i.e. the steady state 
distribution given the concentration Thus we get 

La;{x,x) = Li{x,x,Co{x)) (24) 
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and 


P{nf,T\ni,0)(x J I?a;exp-|— J dtLx{x,x) 


(25) 


Although we still leave the factor P “ in the Laplace aysmptotics (231 in our manipulation and then take the singular 
perturbation analysis, it is not difhcult to establish the final result in a rigorous way. This result tells us that 
when a > 1 the LDT for the slow variable x is only determined by the effective synthesis rate 
degradation rate 7 . 

(ii). Case 3: 0 < a < 1. The translation process is much faster than the switching process. 


In this case, we rewrite Eq. (191 as 


P(u/,T|ni,0) oc J I?xZ3c exp I — ^ y dt ^ ^_^ Li{x,x,c) + L 2 {x,c)^^ . 


Taking Laplace asymptotics with respect to c-integral, we get 


Lx{x, x) = inf ^-^Li(x, x, c) + L 2 {x, c)^ . 


(26) 


(27) 


Since Li > 0 and e°‘~^ ^ 1, the term e““^Li(a:, i, c) dominates. We can perform similar approach to derive as in 
the previous case. In general, we assume that Li{x,x,c) achieves its minimum Li(x,x,c) = 0 at c = for a given 
X. In our case, c^, satisfies the mean field ODE by large volume limit: 


N 


By (271 and (28), we have 


and 


* = 

j=i 


Lx:{x,x)= inf L 2 {x,Cx) 

{cx'x='^^_^ kj{cx)j—'yx} 


P(n/, rjui, 0 ) oc J i7a:exp|— -J dtLx{x,x) 
We want to remark here that from (30) we will expect to get the LDT of the type 

lim e“lnProb(X S S) = — inf [ Lj;(x,x)dt 

€.—^ 0 + X^ISJq 


(28) 


(29) 


(30) 


(31) 


where S is a Borel set in D[0, r] space (functions on [OyH are right continuous with left limits) and X. is the sample 
path of the original jump process. The scaling e“ in ( |31[ ) is essential to reveal the nontrivial behavior of x. Other 
choices of the exponent do not give the correct limit which we are interested in for x. 

In the considered case 0 < a < 1, the protein synthesis are much faster than the genetic switching. With this 
condition, if we neglect the copy-number fluctuation of the protein, we get a reduced stochastic hybrid system: 


N 




(32) 


where /{ 5 (t)=j} is an indicator function and ^{t) represents the DNA occupation state. In [3] and [ 8 ], the authors 
established the LDTs for variable x as e ^ 0-1- for the system (32), which is the same as what we derived in Eq. 
(29). But we should emphasize that this coincidence is not obvious a priori, our result supports the validity of the 
procedure by taking the repeated limit for two-scale processes in some sense. 

(iii). Case 2: a = 1. The switching rates are comparable to the translation rates. 

When a = 1, we have 


P(n/,T|ni, 0 ) = J DxDcexp^^-- dtLi{x,x,c) + J dtL 2 {x,c) 


(33) 
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In this case, we have the LDT Lagrangian for variable x: 

Lx{x,x) = inf {Li(a:,x,c) + L 2 {x,c)} . 


(34) 


Since in most cases there is no closed form for Li, thus we do not expect to get the closed form of accordingly. 
This hinders the applicability of the obtained theory. It is more convenient to study the conjugate Hamiltonian of L^: 

Hx{x,p) = sup {p/3 - Lx{x,l3)} 

P 

= sup - inf {Li(x, /3, c) + ^ 2 ( 2 ;, c)}| 

= sup sup {p/3 - Li{x,P,c) - L 2 {x, c)} 

P c 

= sup sup {p/3 — Li{x, 13, c) — L 2 {x, c)} 

c /3 

= snp^Hiix,p,c) - L 2 ix,c)Y (35) 


As we will show, the dual Hamiltonian may have explicit expression and it is convex with respect to the momentum 
variable p. This property makes it competitive for the nu merica l algorithms for solving static Hamilton-Jacobi equation 
through the geometric minimum action method (gMAM)(SE3. 


IV. APPLICATION TO THE TWO-STATE MODEL 


Using the two-state model 0 as an example, we will give the detailed LDT results for different a, and show the 
mean field ODE and the chemical Langevin approximation for variable x. Moreover, we will solve the static Hamilton- 
Jacobi equation for the quasi-potential <I>(x) in different situations. At first, we take the same rescaling (13) for the 
variables and parameters. We again consider three different cases: (i) a > I, (ii) 0 < a < I and (in) a = 1. 

(i). Case 1: a > I. 

In this case, the ergodic limit of DNA occupation probability is 


Co(x) = 


fix) 


gix) 


for given x. By Eq. (23), we have 


Lx{x, x) = sup Ipx — 

p 


fix) + gix) fix)+g{x) 

kifix) + hgix) 


fix)+g{x) 


(exp(p) - 1) -b 7 x(exp(-p) - 1) 


and the dual Hamiltonian 


From the result 


we get the mean field ODE 


Furthermore, the fact 


Hxix,p) = ^ [exp(p) - I] -b 7 x[exp(-p) - 1]. 

fix) + gix) 


dHx _ kifjx) + kpgjx) 
dp p=o Jix)+g{x) 

dx dHx kifix) + kogix) 


IX, 


dt dp p=o fix) + gix) 

d'^Hx _ kifjx) -b kogjx) 
dp'^ p=o /(x)-bp(x) 


"fX. 


■ 7X 


(36) 


(37) 


(38) 


(39) 


(40) 
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shows the following chemical Langevin approximation holds 


cfc = (hihllMM _ ^ 

V fix)+Hx) j 


lkif{x) + koKx), ,— , \ 

/M + 5W . 


(41) 


where Wi and W 2 are independent standard Brownian motions. 

From classical variational analysi^^, it can be shown that the quasi-potential defined through 


$(a:;xo) = inf inf / Lx(x,x)dt 

T >0 x(0)=xo,x(t)=x Jq 


in our case satisfies a static Hamilton-Jacobi equation H{x, dx^) = 0, where xq is a stable fixed point. Based on (371, 
we have by some algebra 


dx^ = - log 


(kj(x) + kff(x))/(f(x) + g{x)) 


7X 


(42) 


This result is consistent with the quasi-potential derived in m, where the authors neglect the fluctuation of genetic 
switching and get the result by WKB ansatz. But of course, there is no hope to get the explicit formula of $ when 
the dimension of x is bigger than 1. 

(ii). Case 3: 0 < a < 1. 

In this case, we have the Lagrangian 

Lx{x,x) = L 2 {x,Cx) 

= sup |- ^Cig{x){e'^^~'^^ - 1) + C 2 f{x){e‘^^~'^^ “ } 




(43) 


where Cx = ( 01 , 02 ) and Ci = (x — kg + ^x)/(ki — kg), C 2 = 1 — ci by the condition Li(x,x,Cx) = 0. With the 
Legendre-Fenchel transform defined by Hx{x,p) = sup^ (pP — Lx{x,/3)), we get the dual Hamiltonian: 


Hx(x,p) = pPo - Lx(x, l3o), 


(44) 


where /3o = kiSi + k 2 (l — si) — jx and 


Si = X 


S2 


2 ' 27sfT4’ 


S2 = 


p{ki - ko) + fix) - g(x) 


\J f{x)Kx) 

Again, we can obtain the deterministic mean field ODE as 

dx dHx kif(x) -I- kog{x) 


Similarly, we get 


dt dp 


d'^Hx 


P=o f{x)+g{x) 

2f(x)gix) 


7X. 


(f(x) + gix))^ 

and thus the chemical Langevin approximation 


(fci - ko)" 


dx = I MhllMM _ J * + 

V fix) + gix) I 


' ‘2fix)gix) 

ifix)+Hx))^ 


(ki — ko)'^dw. 


(45) 


(46) 


(47) 


We note here that the fluctuation term has the strength ve“ since it originates from the fast genetic switching process, 
and the term yf^dw 2 disappears because it is in order y/e. These are in sharp contrast with the result in (41) and 
Case 2 below which has the 0(\/e) fluctuation. 

























By Eq. (441, solving the Hamilton-Jacobi equation = 0, we have 




/ 


( 48 ) 


kt) — qcc ki — 'yx 

This is consistent with the result in m although we have totally different form of Hamiltonian H. 

(iii). Case 2: a = 1. 

In this case, the genetic switching rate is comparable to the protein synthesis rate. The rigorous LDT result has been 
obtained in m_ Now we formally establish the LDT again through the second-quantization path integral approach. 
By (34) and ( [3^ , we have the dual Hamilton; 


c{x,p) = sup^Hi{x,c,p) - L 2 ix,c)'^ 

= sup I (jiiCi + koC 2 ^ [exp(p) -1]+ ^x[exp{-p) - 1] - - ^/g{x)cl^ 

Since ci -I- C 2 = 1 and Cj > 0, we can obtain the explicit expression of Hx{x,p): 

Hx{x,p) = (kis + Ao(l - s)j [exp(p) - 1 ] + 7 x[exp(-p) - 1] - (^\/ (/(I - s) - ■ 


(49) 


(50) 


where 


s = 


Si 


2/ 


sf -h 4 


Si = 


(fci - ko){eP - 1) -b fix) - gix) 


As before, we get the mean field ODE 

dx dHr. 


dt dp 


'fix)gix) 

kif{x) -b kogix) 


5=0 fix) + gix) 


7X. 


The second order expansion to p 


dp"^ p=o ifix) + 5 ( x ))3 
yields the following chemical Langevin approximation 

fci/ + 


d'^Hx 2fix)gix) ~ ~ 2 , hfix) + hdix) 

(fci -ko) -b ' 


fix) + gix) 


■ 


dx = 


f + 9 


— yx dt -b -s/e 


kif + kpg 
f + 9 


;ih - h)^dwi - y/^d'W2\ 


if+ 9)' 


(51) 


(52) 


(53) 


where f,g are abbreviations of /(x) and gix), and Wi,W 2 are independent standard Brownian motions. 

It is worth discussing the relationship between the Hamiltonian (ISOl) and that obtained by WKB asymptotics. To 


get a Hamiltonian via WKB asymptotics, we follow the procedures in |21j and sketch its outline. We assume that the 
stationary solution of ([^ has the form 


— <i)(x) 


i = 1,2. 


(54) 


P/(x) ~ rj(x)/c(x) exp 

Substituting (54) into (H and collecting leading order terms, we get M(x,p) • (xi(x), r 2 (x))'^ = 0, where 

^ ^P> [ ~g fco(eP-l)+7x(e-P-l)-/ 

Now there is a subtlety to get the correct Hamilton in LDT. If one takes the choice that the iL(x,p) is defined as the 
largest eigenvalue of M(x,p), it can be sh own that it is equivalent to (50). However if one takes another choice that 
it is defined as the determinant of M(x,pP2iJ, then H is not convex in momentum variable p and the equivalence is 
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lost. This resembles the issue of the choice of Hamiltonians for parametrized curve problem in classical mechanics 
(see p. 40 in [5]). 

Let us make some comments on the obtained mean-field limit and Langevin approximations. Recall that there are 
two parts of noise in the original dynamics: one is from the translation process and the other from DNA switching 
process. Our result tells that different diffusion approximations arise according to the magnitudes of residual noise 
in different reaction channels. If a > 1, the dominant part of noise is from the translation, so only the fluctuation 
from translation process survives. If 0 < a < 1, the dominant part is from switching process, so only the fluctuation 
from DNA switching process survives. And if a = 1, both fluctuations from protein translation and genetic switching 
contribute. Similar situation occurs in the LDT analysis, where the singular perturbation is performed for variational 
minimizations. The obtained results show the validity of the procedure by taking the limit for faster process at first 
and then performing the corresponding analysis for slower process. Although we only consider the two-state models, 
the essential structure and results hold for general cases. It is a natural extension of the classical large volume limit 
for chemical reaction processes. We summarize our discussions for the three regimes in Table |Tj 

TABLE I: Comparison of LDTs, mean-field limits and Langevin approximations for three regimes 



LDT Hamiltonian 

Deterministic Drift 

Noise in Langevin Approximation 

Case 1 : (a > 1) 

H^{x,p) in (37l 

kif{x) + kog{x) 

— 

f{x)+g{x) 

A'l 

f + g J 

Case 2: (a = 1) 

Hx{x,p) in (50l 


H- — (fci — koYdwi — •J^dw2 1 

{f + gY ) 

Case 3: (a < 1) 

H^{x,p) in Q 




V. APPLICATION TO THE SINGLE-MOLECULE MICHAELIS-MENTEN KINETICS 


Our approach and observation have interesting implications on the single-molecule Michaelis-Menten systeirfi^, in 
which a substrate S binds reversibly with an enzyme E to form an enzyme-substrate complex ES that decomposes to 
form a product P. The reaction schemes can be schematically shown as 


fci ^2 

E-hS .. .> ES -► E-bP, 

fc-i 


(55) 


In case of single-molecule enzyme set-up, the reaction system (551 falls in the framework considered in this paper. As 


in [HI, we assume that the substrate is abundant enough and there is essentially no depletion of substrate by a single 
enzyme molecule. That is, we assume the concentration of substrate is a constant, which will be denoted as [S']. It is 
well-known that the rate of product formulation v has the following form in the quasi-steady state approximatiorP^ 


k2[S] 

[S] + kM 


(56) 


where kM = {k_i + k 2 )/ki. In [H], the statistics of enzymatic turn-over time and dynamical disorder are considered. 
Here we are interested in deriving the Langevin approximations of the Michaelis-Menten system in different regimes. 

In [H], k-i ranges from to 2000s“^, fci is usually taken as ^2 = 250s“^, and [S] ranges from the 

order O.OOlmM to O.lmM, where IM = Imol/L. Some specific choices of these parameters include 

• Case 1: fc_i = 20005"\ fci [S'] = lO^M-^s-^ x O.SOmM = 3000s-\fc2 = 250s-i; 

• Case 2: fc_i = 200s-i, A:i[S] = x 0.02mM = 200$-'^, k 2 = 2505-^; 

• Case 3: fc_i = 50s-\ fci[S] = lO^M-is"! x O.OOSmM = 50s-i, fcz = 250s-i. 


10 





























The degradation rate constant 7 is not essential and we assume it is 0(ls ^). The above choices underlie the rationale 
to study different regimes in previous sections since we can make the assumption 


ki[S],k-i ^ a>0 


(57) 


if we define e = 1/250. 


The chemical master equation of the system (551 can be written as 

/ 


dtP{n,t) = 


0 fc. 


(E„ - l)P(n,t)+7(E„ - l)[nP(n,t)] + 


0 0 


^-ki[S] /c-i+fca 
ki[S] -(fc.i+fe) 


P{n,t). (58) 


Denote ci the occupation probability of the free enzyme molecule state E and C 2 the probability of the complex state 
ES. With similar approach in deriving Q, the transition probability can be obtained with Lagrangian L 

- P(n,c,i/3,i(/>), (59) 

where the Hamiltonian 

H{n, c,il3,icj}) = k2C2[exp{i/3) - 1] + 7 ri[exp(-i^) - 1] + C 2 (fc_i + k 2 e'‘^){e"^^ - 1) + ciA:i[S'](e“®‘^^ - 1). (60) 

According to ( |5^ , we make the rescaling 

a: = ne (or Xq, = ne“), ^2 = ^ 26 , ^1 = [<S']e“, fc_i =/c_ie“. (61) 

Then the Cases 1, 2 and 3 correspond to a > 1, a = 1 and 0 < a < 1, respectively. Next let us study the three 
cases separately. The order of discussion will be from easy to difficult, which may be slightly different from previous 
sections. 

(i). Case 2: a = 1. 


With the steepest descent asymptotics as in (191, we have 


P(n/, r|ni, 0 ) cx J DxDcexp J L{x,x,c)dt^ . 


(62) 


The Lagrangian L has the form 

L{x, x,c) = supsupjpx — H{x,c,p, If)} 

p V 

= sup{px - H{x,c,p)}, 

p 

where H{x,c,p) = inf^{P(x, c,p, (^)} and 

H{x, c,p, if) = fc 2 C 2 (e^ - 1) + - 1) + £ 2(^-1 + fc 2 e^)(e‘^ “ 1) + Cifci(e“‘^ - 1). 

In this case, we have the LDT Lagragian for variable x by applying Laplace asymptotics 

La;(x, x) = inf {L{x, x, c)} . 


(63) 

(64) 

(65) 


Following the approaches in deriving (351, we get the conjugate Hamiltonian of L^: 


,(x,p) = sup {p/3 - L^(x,/3)} = sup{p/3-inf{L(x,/?,c)}| 

/3 /3 c J 

= sup sup {p/3 — L(x, (3, c)} = sup sup {p/3 — L(x, (3, c)} 

/3 c c /3 

= sup {P(x, c,p)} = supinf{P(x, c,p, ip)} 

c c ^ 


= k2s{e^ — 1 ) + 7a^(e ^ — 1 ) — y fci(l — s) — v (fc_i + ^26^) 


( 66 ) 

(67) 
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where 


_ £ 5i 

2 2y/sf+A' 


Si = 


ki — (fc_i + £ 2 ) 
^Jki{k_i + fcaeP) 


We have the mean field ODE by local analysis 

dx dHj. 


dt 


dp 


fcifc 


p=o ki + k-i + k 2 


"fX. 


(68) 


This is consistent with the Michaelis-Menten law shown in ( 56 ). To see this, we first note that the reaction rate v 
should be rescaled with e~^ since (68) is for the concentration variable x instead of n. We have 


lim 


k2[S] 


= lim 


kie ^k2e 


k-ih 


l't2 


e-lO+ [S'] + kM e-lOS- 1 -(- k-iE ^ ^ + fc-l + ^2 

Furthermore, the second order expansion of with respect to p 

2km kk 


dp'^ P=o (ki + k-i + k2Y ' ki + k-i + k2 

yields the following chemical Langevin approximation 

kik 


■ jx 


dx = 


k\ + k—i + £2 


— ^x \ dt + y/e 


kik2 




A;i + fc_i+fc2 (fci + A;_i + ^ 2 )^ 


dwi — y/^dw2 ■ 


( 69 ) 


( 70 ) 


Although the above result is quite natural based on our derivations in previous sections, the application to Michaelis- 
Menten system again tells us that the strict correspondence between the drift and diffusion terms in classical large 
volume limit is lost. 

(ii). Case 1: a > 1. 

In this regime, ki[S] and k-i are much larger than fc2- Similar as in Case 2 , we have the transition probability 
density ( 62 ) with Lagrangian ( 63 ), and thus the Hamiltonian (66) for the slow variable x. We get 


^^a;(a:,p) = supinf I fc2C2(e^ - 1)-b7a;(e ^ - 1)-b Ca f -b ^26^ j (e"^ - 1)-b Ci ^^(e - 1) 


= sup { k2C2{eP - 1) -b -fx{e p - 1 ) - f 


fc_i -b ^ — V ci^i 


( 71 ) 


As e —>■ 0 -b, the singular perturbation analysis suggests the term involving l/e“ ^ to be 0 , which gives ca = fci/(fc_i -b 
ki). We obtain 


H^{x,p) 


r m - 1) + - !)• 

ki -b /c_i 


( 72 ) 


The mean field limit 


dx 

dt 


kik2 

fci -b fc-i 


— 7a;. 


Its consistency with the Michaelis-Menten law is straightforward by checking 


lim e 
£->■ 0 - 1 - 


k2[S] 

[S'] -b kM 


lim e-^r-- 


fcie “fcae ^ 
-b k-ie~°‘ + 


fcifca 

ki -b k—i 


(73) 


for a > 1. 
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In this regime the chemical Langevin approximation takes the form 


dx = 


kik2 


— \ dt + 


~ — dwi — y/^dw^ 

ki + fc_i j 


k\ + k—i 

and we formally recover the correspondence between the drift and diffusion terms in this case. 

(iii). Case 3: 0 < a < 1. 

In this regime, we need to select the scaling parameter for n as e“, namely to define Xa = ne°‘. We have 


(74) 


P{nf,T\ni,0)o(. J DxaDcexp^—-J L{xa,Xa,c)dt 


(75) 


where 


L{Xa,Xa,c) = SUPpSUp,^ ^pXa “ ^I3^C2(e^' - 1) -7X0(6 - 1) 

-62 (k-1 + ] (6^ - 1 ) - Cifci(6"^ “ 1 ) r • 


(76) 


The Hamiltonian corresponding to variable Xa has the form 

Hx{xa,p) = sup inf I -^^C 2 {eP - 1 ) + 7 X 0 ( 6 “^ - 1) + 62 [l-i + j - 1 ) + ciAi(e“‘^ - 1) 


= sup < -^C2{eP - 1 ) + 7 X 0(6 P - 1 ) - 


A 


k2eP 


C 2 *-i + - V 


As 6 —>■ 0+, the singular perturbation analysis gives 62 = °‘kieP/k 2 , and the final Hamiltonian 

Hx{xa,p) = ki{e^ - 1) + 7a;a(6“^ - 1). 

We get the mean field ODE 

dxa r 

— = ki- 7 x 0 , 

The consistency with the Michaelis-Menten law can be verified by checking 


lim 


k2[S] 


kie “^26 ^ 


= lim 


£->■0+ [<S'] +/cm £->-0+ kiC °‘ + k-ie “ + ^26 ^ 

for 0 < a < 1. In this regime we have the chemical Langevin approximation as 


= ki 


dXa = (j^l — 7 X 0 , j dt + y/e^ — y/jXadW 2 ^ . 


(77) 


(78) 


(79) 


Again the formal correspondence between the drift and diffusion terms is recovered but with special concentration 
variable definition and rescaling. It is instructive to compare (79) and (47) in the case 0 < a < 1. We have an 


additional term y/jXadw 2 in (79) because of the utilized scaling Xa = ne°‘ instead of x = ne in the two-state model. 
This reveals the difference between the single-molecule Michaelis-Menten and the two-state genetic switching model. 
Indeed in this regime, the reactions can be simplified to 

ki 7 
S--^0 


since the production rate is limited by the rate of formation of complex ES. 

We summarize our findings for the single-molecule Michaelis-Menten kinetics in three regimes in Table [ill 


13 





















TABLE II: The LDTs, mean-field limits and Langevin approximations for single-molecule Michaelis-Menten kinetics 



LDT Hamiltonian 

Deterministic Drift 

Noise in Langevin Approximation 

Case 1: (a > 1) 

Hj:{x,p) in (72 1 

k\k2 

lAI 1 / 

■J 

1 

1 

+ 

1-^ 

ki + k-i J 

Case 2: (a = 1) 

Hj:{x,p) in (67l 

k\k2 

\^(J k\k2 2k^k2 ^ 

ki k—i /t2 

\\j ki + k-i + k2 (fci -b k-i + k2)^ 

Case 3: (q < 1) 

Ha:{x,p) in (77l 

fei - 'yXa 



VI. DISCUSSIONS AND CONCLUSION 


The methods and LDT results we proposed in this paper are not limited to the two-state model, single-molecule 
Michaelis-Menten and single kind of product case. It is indeed general for a class of two-scale kinetic systems. To 
show this, let us consider the following extension as shown in Fig. 


State7 


f{ni,n2] 


State2 



Product7 {ni) 


g('ni,n2) 


^22 


Product2 (n2) 




FIG. 2: Schematics of a two-scale kinetic model with two kinds of products. 


We assume similar scaling as considered in (12): 

1 
e 


h- ■ 


a>0, {i,j = 1 , 2 ). 


Define Xj = rije, kij = kijC for i,j = 1,2 and f(xi,X 2 ) = fini,n 2 )e°', g{xi,X 2 ) = g(ni,n 2 )e°‘■ Performing the same 
approach as in Sec. |Tlj we get the transition probability 


P(n/,r|ni,0) cx 


DxDc exp 



dtLi{x, X, c) 


1 

e“ 


dtL2{x, 



Here the Lagrangian 


Li{x,x,c) = snp{p ■ X - Hi{x,c,p)}, L 2 {x,c) = sup{-H 2 {x,c,(p)}, 

P V 


where x = {xi,X 2 ), c = (ci,C 2 ), p = (pi,P 2 ), = { 9 ^ 1 , V 2 ) and 

Hi{x,c,p) = ( 7 iici -b fc2iC2)(ePi - 1) +'yixi[e~^^ - 1] + (fci2Ci -b fc22C2)(e^" - 1) + 722:2(6“^'' - 1), 
H2{x, c, if) = cig{x){e‘^^~‘^^ - 1 ) + C2f{x){e'^^~'^^ - 1 ). 


All of the analysis performed for the two-state model can be applied here to obtain the LDTs for variable x = {xi,X 2 ) 
with different a. 

One can also employ the WKB ansatz Pj{x) ^ exp(—for the stationary distribution of the stochastic 
hybrid system (321, where x = ne and j is a state of DNA. In the asymptotics, one gets a static Hamilton-Jacobi 
equation for the quasi-potential and it turns out does not depend on the specific choice of L However if not 
handled appropriately, the WKB approximation may lead to totally different forms of HamiltoniarPSl as mentioned in 
the end of Section |IV[ This non-uniqueness is due to the lack of variational selection in LDT, which gives a unique 
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Hamiltonian dual to the obtained Lagrangian in rate functional. And this Hamiltonian has the superiority that it 
is convex with respect to the momentum variable as the by-product of Legendre-Fenchel transform and LDT. This 
property is important for the nice behavior of numerical discretization. 

In this paper, we assume the switching rates between different DNA states are in order e~°‘. It is not necessary 
and could be more general. As long as the switching rates between different DNA states are in 0(A(e)), and the cases 
lim£_j.o-i- eA(e) = oo, 0(1) and 0 are considered, we will get similar results. Especially, the readers may easily verify 
that if we assume A(e) = Ke~^, then the two-scale LDT Lagrangian with e-scaling in front of — InP has the form 

Lx{x, x) = inf {Li{x, x, c) -|- KL 2 {x, c)} , (80) 

C 

and the LDT Lagrangian with A(e)-scaling in front of — InP has the form 

Lx{x, x) = inf {K~^Li{x, x, c) + L 2 (x, c)} . (81) 

When K goes to 0,1 or oo, the appropriate choices of scaling recover the desired results shown in the paper. 

In conclusion, we established the two-scale LDTs for a class of chemical reaction kinetics through the second- 
quantization path integral formulation. Although not rigorous, we showed that this formal approach is very effective 
and transparent to understand the two-scale LDTs associated with different reaction channels. This provides essential 
insights to rigorously prove the corresponding LDTs, which is our ongoing research. We discussed its implication 
on single-molecule Michaelis-Menten kinetics as well. The proposed framework and results also shed lights on the 
understanding of general multi-scale systems including diffusion processes. It will be interesting to investigate the 
application of two-scale LDTs to other systems. 
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